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Abstract; The motion of a satellite in orbit, subject to atmospheric force 
and the motion of a reentry vehicle are governed by the same forces, namely, 
gravitational and aerodynamic. This suggests the derivation of a uniform 
set of equations applicable to both cases. 

For the case of satellit** motion, by a proper transformation and by 
the method of averaging, a technique appropriate for long duration flight, the 
classical nonlinear differential equation describing the contraction of the 
major axis is derived. WTiile previous authors, and in particular King-Hele, 
integrated this equation using various heuristic methods, the present authors 
present a rigorous analytic solution, with a high degree of accuracy, using 
Poincarfe' s method of small parameters. Next, using Lagrange's expansion, 
the major axis is expressed explicitly as a fxinction of the eccentricity. The 
solution is uniformly valid for moderate and small eccentricities. This is a 
major achievement due to the discovery of a certain recurrence formula which 
facilitates the long and tedious analytic process. For highly eccentric orbits, 
the asymptotic equation is derived directly from the general equation. To 
obtain the same equation King-Hele must use an entirely different method. 
Again, while King-Hele only succeeded in obtaining an approximate solution 
to this case using a heuristic method of integration, the exact solution to the 
asymptotic equation has been obtained by the present authors. Nximerical, 
solutions have been generated to display the accuracy of the analytic theory. 
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INTRODUCTION 


The theory of satellite orbits in the presence of an atmosphere was 
developed during the late fifties with the launching of the first artificial 
satellite. With increasing knowledge of planetary atmospheres, especially 
the atmosphere of the Earth, the theory has now reached a very high degree 
of accuracy. For a first estimation of the lifetime of a satellite and for a 
correlation between the semi-major axis and the eccentricity of the orbit 
while it undergoes a contraction due to the perturbing effect of atmospheric 
drag, analytic theory is adequate. The classical theory was presented in 
a monograph written by King-Hele (1964), who is among the authors who 
have contributed the most to the development of analytical solutions. 

There are two reasons for presenting this new study of a well 
established and analyzed subject. 

The first reason concerns the approach to this problem. In the 
early days, development of the theory of flight inside an atmosphere was 
conducted from two different approaches. On the one hand, researchers 
analyzed the small perturbations of satellite orbits at very high altitudes. 

The mathematical tools are perturbation theories in celestial mechanics 
based on Lagrange's equations for the variations of orbital elements. TI.c 
space object, usually referred to as a satellite, is not intended for recovery. 
The main subjects of concern are its life expectancy and, while in orbit, 
the slow variations of its orbital elements. On the other hand, engineers 
and scientists who were concerned with the safe recovery of an entry 
vehicle concentrated their effort in the study of the deceleration and heating 
during entry. Here the elements of prime consideration are the position 
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and velocity of the vehicle, both varying rapidly. The smooth behavior 
of near Keplerian orbit is no longer valid and strong physical assvimptions 
were made to such an extent that, although describing the same phenomjnon, 
namely flight of an object inside a planetary atmosphere, the equatiois 
became totally different. The gap got wider, as the two theories became 
more and more sophisticated, so that now the two groups, one consisting 
mostly of mathematicians, and one consisting mostly of space dynamicists 
seldomly reference the other group's work. With the objective of providing 
a unified theory for flight inside a planetary atmosphere, w’e have formulated 
a set of xmiversal, exact equations. These equations have been successfully 
applied to the study of planetary entry of a space vehicle (Vinh et al. , 1977) 
and even to optimization of such an entry (Vinh et al. , 1975). In this paper 
we shall present the necessary transformation such that the equations can 
be used for analyzing the slow variations of the orbital elements while the 
vehicle is still in near vacuum. This successful wedding is necessary since 
the future space vehicle is designed to stay for an extended period in orbit 
as a satellite, and also to be recovered safely after a fiery entry followed 
by a glide, an approach and a landing on an airfield. 

The second reason concerns an improvement of the existing theory. 
Often, because of mathematical difficulties in the analytic integration of 
the equations of motion, the types of solutions are artificially 
classified. One can easily understand the classification of orbits into 
hyperbolic and elliptic because the natures of the orbits are different, as 
reflected by the Keplerian equations with e > 1 and e < 1 . However, 
when it comes to different orbital phases with e very small and with e not 



80 small, then it is clear that the classification is purely for easing the 
integration and usually the regions of validity of different solutions are 
at best defined empirically. Our eHort in going over the classical theory 
is to remove, whenever possible, such an ambiguity. 

FORCES ON A SATELLITE IN ORBIT 


The satellite and the planet are assiimed to be in two-body relative 

motion. For a spherical planet, the gravitational force is an inverse 

square force of attraction with acceleration 

g(r) = ^ • (1) 

r 

The atmospheric force is in the form of drag, D. , acting in a direction 
opposite to the velocity V relative to the ambient air 
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is the drag coefficient measured using a reference area S and the 

atmospheric density p . We shall use a strictly exponential law 

Mr -r) 

p = p e (3) 

^o 

where p is a constant 

P = ^ (4) 

H is the scale height, and subscript p^ denotes the initial periapsis condition. 

The equations of motion are written with respect to an inertial frame 
with the origin at the center of the planet. Let V be the absolute velocity 
of the satellite. 
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^ \ ^ \ ( 5 ) 

where is the velocity at the point M, ol the ambient air relative to the 
planet center (Fig. 1). If w is the angular velocity of the rotating 
atmosphere, then 

V = r w cos (j> (6) 

where ^ Is the latitude of the point M. 

Let be the angle between and V. Then by squaring Eq. (5) 

+ V ^ - 2 VV cos 4 j' . (7) 

A e e 

The vector V is in the local horizontal plane. Also, near the periapsis 
© 

where the aerodynamic drag is most effective, the satellite travels in a 

nearly horizontal direction and hence the angle y between the velocity V 

and the horizontal plane is small. Then, the angle between and 

V is nearly equal to the angle 4< between V and the projection V of V 

© rl 

on the horizontal plane. This angle 4^, called the heading, it: related to 
the latitude 4> and the inclination i of the orbital plane by the well* known 
relation 

cos cos 4> = cos i . (8) 

Therefore, we have approximately 

V cos O' a V cos 0 = rwcoS(j> cos 0 = rwcos i. (9) 

© © ’ 

Upon substituting Eq.i. (6) and (9) into Eq. (7) we have 

2 ,,2 ,, rw .2 2 2, 2 , 2 

= V (1- -^-cosi) + r w (cos <|> - cos i) . (10) 
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The rotation o£ the atmosphere is generally slow so that the term 
can be neglected. For the small term rw/V, it is appropriate to use an 

average value. King-Hele suggested using the value at perigee r /V 

P P 
o o 

to replace r/ V. Finally i, which usually varies by less than 0. 3° 

during a satellite's life, may be taken equal to its initial value i . Then 

o 

we have King-Hele's expression (King-Hele, 1964) 

( 11 ) 

where the average constant value f is 

r w 

P 2 

f = ( 1 - COS i^ ) . (12) 

Po 

Thus, in terms of the absolute speed, the drag force is 

acting opposite to the direction of the velocity of the satellite relative 
to the ambient air. 


THE EQUATIONS OF MOTION 

For the flight of an aerodynamic vehicle with a lift coefficient 
and a drag coefficient Cj^ , it is customary to use the equations of 
motion with the notation in Fig. 1. 


dt 


V sin Y 


dt 


V cos V cos ill 
r cos (|> 


d6 

dt 


V cos Y sin iji 

r 
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dV 

dt 


= 

dt 


V ^ = 

dt 


pSCj^V^ 


2m 


p SC^ V cos (T 
2m 


g sin Y 


- ( g ■ — ) cos Y 


( 14 ) 


p SC^V sin cr y2 


cos YC 08 4> tan 4> 


2m cos Y 1 

where the bank angle <r is defined as the angle between the local vertical 
plane containing the velocity and the plane containing the velocity and the 
aerodynamic force. 

Using the dimensionless variables 


u = 


Z = 


V cos Y 


gr 


pSC 

2m 


^VT 


(15) 


and a dimensionless independent variable 

‘ V 
o 


r V 

j ( “ ) cos Y dt 


(16) 


we have derived the exact, universal equations for entry trajectories 
into a planetary atmosphere assumed to be at rest (Vinh et al. , 197 5) 


d Z 
ds 

= -P r / - ^ 

Pp 

d£ 1 

dr " 2p r 


Z tan Y 

du 

ds 

Zu 

COS Y 


COS O' tan y + 

sin Y 
2^ Z 

lx 

. z r 


cos Y /I 
+ (1 ■ 

2 

cos Y , 

d s 

COS Y L 

— CO.. 

u 
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coa 4* 


ds 

cos ^ 

(17) 

ds 

s sin 4< 


dv|i 

ds 

Z , cos^Ycos vli tan«t> 

" > • 

cos Y D ^ 


We shall use the necessary transformation to obtain the equations 

for satellite motion inside a rotating atmosphere. 


First, 

still with the notion of an atmosphere at rest, we use 

the 

transfo rmation 


cos <|> 

cos = cos i 


cos ^ 

sin = sini cos » 

(18) 

cos a 

- cos cos ( 6 > Q) 


to transform the last three equations of (17) into 


dg 

ds 

= 1 . Z , ,i„„ 

tan i cos Y D 


dQ 

ds 

Z sin a . , 

= iC ^ ) sin 0- 

sin i cos“ y D 

(19) 

d i 

Vbt Z cos a , , 

= 2 ( ^ ) sin (T 

cos“ Y 


d s 



From Fig. 1 we notice that i is the inclination and JJ the longitude o£ the 
ascending node of the osculating plane. The angle a is the angle betw’een 
the ascending node and the position vector. The angle j- . as stated 
earlier, is the angle between the vertical plane passing through the 



velocity, the ( r , V) plane, and the plane containing the aerodynamic 
force and the velocity, the (A, V) plane (Fig> 2). For satellite motion, 
we have a simplification, and at the same time, a complication. The 
simplification is that there is no lift. The complication is that the drag 
force is modified by the factor f, Eq. (12), and it is directed opposite to 
the velocity , not the absolote velocity V. 

Figure 2 is the aerodynamic force diagram used in deriving the 
eq\iations of motion (14), to which we have added the velocity with 
respect to the ambient air, and the drag force , opposite in direction 
to V . In the present situation, we remove the lift force L and replace 
the vector drag D by D^. This force can be decomposed into one 
component in the orbital plane and one component normal to the orbital 
plane. Since V is small, V is nearly aligned to ^ and the drag 
component in the orbital plane can be considered as directly opposite 
to V, with magnitude D. as given by Eq. (13). To have the component 

A 

of 5^ orthogonal to the orbital plane we find the projection of 


1 2 A 

2 P D V, 


( 20 ) 


By the vector relation (5), since V is in the orbital plane, and since 

V makes an angle 4i with the orbital plane, the projection of V on 
0 A> 

the normal to the orbital plane is the same as the projection of which 
has magnitude 


V sin ^ = rwcos 4) sin4« s 


- 


rw sini cos a 


( 21 ) 


Hence, the vector has magnitude 
1 

= TP Sf C_ rw sin i cos a tt“ 
N 2 D V. 

A 


By relation (11), we write it as 
pSfCj^V 


( 22 ) 


N 


2 £ 


TTT 


Tw sin i cos a 


(23) 


and its direction is opposite to the vector L sin <r in Fig. 2. 

The end result of the analysis is that, in the Eqs. (17) and (19), 
we replace Cj^ by the modified drag coefficient f Cj^ , we delete the 
component Cj^ cos<r and replace the component Cj^ siny by 

1/ 2 rw 

sin 0 - = -f Cjj ( -^ ) sin i cos a . (24) 

Finally, the variable Z, the modified Chapman Z function defined 
by Eq. (15), is most effective in analyzing the entry phase of the vehicle. 
While the vehicle is still in orbit, we use it in the form 


- P(rp-r) 

^ Z s (t?-) e ° 


o r 


(25) 


where the dimensionless constant Z is 

o 

p Sf C_r 

— P« ° 

o o 

z = ^ 

o ■» 

2m 


(26) 


We can now rewrite the Eqs. (17) and (19), introducing the equation for 
r/ r to replace the equation for Z 


_d 
d s 


( 



P 


o 



tan Y 
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2Zu "■» 

o i ^ , O 

- -7STT ‘T-’* 

Po 


s 1 - 


r wZ P(r -r) 

p o 5/2 . . p 

, . o , r , cos 1 sin a coj a o 

* 1 + (- — ) — ni 

p^ u COSY 


r w2 5/2 P(r -r) 

p o . ^ p 

*^o , *■ I sino cos a o 

7==- ’ "172 ' 

ypf/ Tp U CO> V 

r wT , p(r -r) 

p o 5/2 . . 2 ^ ' p 

o , *■ . sin i cos or *^o 

^ 1/2 ® 

^pf/ Tp p^ u cos Y 


THE VARIATIONAL EQUATIONS 

The equations (2 7) are the bridge between satelHce theory and entry 
theory. As a matter of fact, they can be used to follow the motion of a 
vehicle subject to gravitational force and drag force of a rotating planet 
for its entire life in orbit until the end of its entry and contact with the 
planetary surface. The accuracy depends on the readjusmient, for each 
layer of the atmosphere, of the "constant value" p. The equations are 
most useful for analyzing the last few revolutions and the entry phase. 

The variables a, S2 and i which are orbital elements are related to the 
entry elements o , 9 and o through the relations (18). On the other hand, 
the variables r, u and y which are the entry variables can be transformed 
into the orbital elements through explicit relations. 

Consider the osculating orbit, that is, the orbit the vehicle would 
follow if at any time the drag force suddenly vanished. Putting s 0 
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in Eqs. (27), we have 


d 

ds 




tan Y 


du 

d 3 

ix 

d 3 

da 

d $ 

dn 

s 

d i 
d s 


- u tan Y 

2 

1 . S21-J: 

u 


I 

0 

0 


The integration ie sinnple and we have the general solution 

2 


cos Y ■ 


u 


2u - C 


1 


u 


u 


1+ - Cj cos(s-C^) 


= a + C, 


n = c. 


(28) 


(29) 


We see that s is equivalent to a and actually, we have only 5 constants 
of integration. The last constant of integration is obtained through the 
time equation, Eq. (16). In the first three equations, we evaluate the 
constants of integration bv taking the origin of time at the time of passage 
through the periapsis. \Vc have 
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( 30 ) 


2 u 

CO* Y = J~ 

2u . (1 . e ) 

u = 1 + e cca { o - w ) 

r = * 

I + e CO8(0r -U' ) 

These three eqviation* provide the link between the entry variables 
r, u and y SLtid the semi-major axis a, the eccentricity e and the argument 
of the periapsis w , which are the orbital elements used in the theory of 
orbits. 


During the phase in orbit, Z is small and the orbital elements 

o 

vary slowly. By taking the derivati\'es of Eqs. (30), considering a, e and 
u as varying quantities and using the Eqs. (27) for the derivatives of r, 
u and Y we have the variational equations for a, e and « . 


First, for e, we have 


22 u 
o 


e cos Y 


?(r -r) 


We present the equation in this form to show an intere -ting behavior of the 
eccentricity of the csculating orbit. It is a general belief that the eccentricity 
decreases continuou.«^ly under the action of atmospheiic drag. This, 
however, is only the secular effect. During each revolution, the flight 
path angle passes through a maximum and a minimum as seen by the third 
equation of (27), and by Eq. (31) it is seen that, at the same time, the 
eccentricity passes through a minimum and a maximum respectively. 

Next, we shall use the more familiar eccentric anomaly E to 
replace s as an independent variable in the variaticnal equations. The 
following relations can be easily derived. 
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cos Y - 


n - e") 

1 4- a cos (or - u ) ” 

(1 

1 > ecos £ 

( 1 - e^) 

(1 - e cos E)(l -t-e cos E) 


= 1 > e cos E 


^ l-e*- 
1 - e cos E 


Hence, changing the independent variable from s to E, the equation for 
e has the form 

1/2 

/I 2 a . _ , l + ecosE , *^o 

dE ' - )(— < ri ' c ^ E -* ' • *”> 


Similarly, we have for the semi-major axis 


da — a 

~ = - 2 Z 
dE or 


( 1 + e cos E)' 
(1 - ecos E’/ 


P(r -r) 
^o 


The last three equations of (27) become 


o I a 


d E d £ 1 - e cos E 


— I sin E(l-e cos E) 


P(r -r) 

1/2 

X (1 + ecosE) e 


r w Z 

Po ° 


|if/ r Vl - Pq 

Pq 


5/2 5/2 1/2 

j-“ j (1-ecosE) (1+ecosE) 


X sin a cos a e 


p(r -r) 
^o 


r w Z 


di_ ^o^ ,a\ _] 

■Tr = “ ~~ I; — I (1-ecosE) (1+ecosE) 

Ifif/r 11-e *^o 

Po P(r -r) 

• • 2 
Sin 1 cos a e 
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THE CONTRACTION OF ORBITS 


The Average Equation 

Consider the variations in e and a, Eqs. (33) and (34). Under the 
dissipative effect of the drag, the naajor axis decreases while the eccentri- 
city, although having an oscillatory behavior, also decreases secularly 
with the time. 

For the radial distance we have, 
r = a (1 - e cos E) 


(36) 


r = a (1 - e ) 

p o o 

* 


We write the exponential function in the equations 


[P (r -r)] = exp [ p (a -a-a e ) + 3 a e cos E ] . (37) 

Po o o o 

Along each revolution a is nearly constant while the varying quantity 
p aecosE provides the fluctuation in atmospheric density. This leads to 
a natural choice of the dimensionless variable 


X = p a e 


to replace the eccentricity. The equation for e is replaced by 

1/2 


dx 

dE 


^^o^* , 1 +e cos E , 

— (e + cos E)( - — rrrr^) 


Mr -r) 


1 - e cos E 


(38) 


(39) 


The new variable x behaves like the eccentricity e; that is, during 
each revolution x passes through stationary values when cos E = -e. On 
the average, however, x decreases with the time. Since the decaying process 
is slow we can use the averaging technique (Bogoliubov and Mitropolsky, 1961), 
applied to the rig’nt h. T’ld side of the Eqs. (34) and (39) for a and x. 
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For the equation for a, we have 


= - 2Z exp[(3(a -a-a e f 

dE or *^'■^'0 oo'-'Ett*' 


2ff „ ^,3/2 

(1+e cos E) 


o (1-e cos E) 


1/2 


X exp (x cos E) dE 


( 40 ) 


For small eccentricity, the integrand can be expanded in a power series 


in e and upon integrating, we have 
2Z a^ 


_ exp [ p (a -a-a e)] |l + 2eI +-je^(I +I ) 
dE r *^‘•'^'0 0 0^)0 1 4 o2 

P. ( 


3 4 


(41) 


where I (x) is the Bessel function of the first kind and of imaginary 


argument, of order n 


1 

I^(x) = J cosnE exp (x cos E) dE 


(42) 


Similarly, the average equation for x is 


dx 

dE 


- 2Z 8a 

O' 


exp [ 0 (a -a-a e ) 1 
^ ‘ o o o •' 


Ij 4 ^ e (31^41^) 


" T'‘”i*V^T6 


4 ,^3 (78 1, 4 3U3 4 3I3) 4 O(e') 


(4 3) 


Let 


Z = — 
a 

o 


(44) 


be the dimensionless semi-major axis. By dividing the Eq. (41) by Eq. (43) 
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and expanding the ratio in a power series in e we have 

^ + r 2y^(3y^+y,)^-29y_^-2y^.y^y J 

3 


o 'o'Z' 8 


'o'-''o'^2' 


'o "^2 ^0^3' 


+ [-32+113y^^38y^y2-y^y^+2y^+6y^^3+2y^y2y3-2y^(3y^+y2)^] 




+2y (11+yJ +8y (3y ^y )(7y +8y +y ) + 1 6{3y +y )(19+y J 

o j oObO^e o c s 


^(y^+y^Hll+yj) - y^C" 8+31y^+3yg) - 2(35y^+36y^+y^) ] 


+ O(e^) 


where we have defined the ratios of Bessel functions 

I 


n 


n 


1 


n 1 


(45) 


(46) 


The equation (45) has been derived by King-Hele, to the order e . 

For X > 3 he integrated this equation by using the asymptotic expansions 

of the function y (x). In this case the right hand side of Eq. (45) has a 
n 

very simple form and the semi- major axis a is obtained by quadrature. 
Mathematically, the method of integration is not rigorous since on the 
right hand side of Eq. (45), the eccentricity is a function of a and x by 
the definition (38); hence, the equation is actually a nonlinear equation 
in a. We shall arrange the equation in a form that allows Poincare' s 
method of perturbations to be applied (Poincare^ I960). 

The Bessel functions satisfy the recurrence formula 


I (x) - I ,(x) = — I (x) 

n- 1 n+1 X n 


(47) 


so that any y (x) can be expressed in terms of y (x) and x. For example 
n o 
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* y. 


2 

X 


4 - i y 

2 X o 

X 


-8 48 

— - — ^ -f y +24 
X 3^0 
X 

1+ 2i t 21 . 
i 4 

X X 


X 

12 

T^o- 


192 

“T 


( 48 ) 


Next, the eccentricity e, expressed in tterms of Zand x is 


X 

e = C — 
Z 


(49) 


where 




(50) 


is a small quantity of the order 10 “ or less. Then, we can write Eq. (45) 

dZ ^ 2 X ,, , 2 , ^o , 

:5“ = cy +€ -:5(2-2y + — ) 

dx o Z o X 

1 2 y 2 , , 

^ ^ , a - O o ^ lx 

+ € 5 ( - 8y - ( + 8y + - ) 

,«2 ^o X o X 

a# ^ 


4 x' 
2Z' 
4 

T 


y« y« 

Os.*’ 


X 

y. 


+ € ‘ ^ ( -4 + ZOy " - 10 ~ + 4 -4 - 5-~ + 20— . 1 by * + - ) 
^^3 o X 3 2 X ^ 


4 1 


X 


+ e (32y . 96y ^ + 82 -2- . - . 17 ~ - 24 

,^4 ' X X 2 3 3 

4Z XX X 

3 4 

y y 

+ 4*^ . 16 - 104 — + o4y ^ ) + 0(€ . (51) 

X X 


We see that the true nature of the equation is a nonlinear equation. 

Since t is very small we need not go further with the expansion, and to the 
5 

order c included, the solution of this equation can be considered as the 

4 

exact solution of Eq. (45) truncated to the order e included. We shall 
use Poincare's method of small parameters for the integration of Eq. (51). 
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Integration by Poincare' a Method of Small Parameters 

Poincare's method for integration (Poincare^ 1960) of a nonlinear 
differential equation containing a small parameter is a rigorous mathematical 
technique, proven to be convergent for small values of the parameter c • 

It has been used extensively in analytic work in celestial mechanics 
(Moulton, 1920). 

We assume a solution for Z of the form 

Z s Z^+eZj +e^Z2 + €^Z2 + €‘*Z^ + €^Zg + ... (52) 

Upon substituting into Eq. (51) and equating coefficients of like powers in 
C , we have the equations for the 


dZ 


1 


dx 
dZ 
"dx 

dx 

dx 

1ft 

dx 


° , 2 

2 1 3 

' — J- (2 + ~ - 2y ) + J ^ "7" ^ 

o o 

2, 2 7 2 T 

Z o Z o 

° 2 ° 

+ -JL_.(- 4 +ij + 20 y 10 — + 4 ^ - 5^+20 1 by '*) 

3 2 o X 3 2 X 'o 

2Z X XX 


dZ, 

dx 


Z ' 7 2 

o Z 


z^’ 

o 


, 7 2 _ 2 
2 , Z, Z, y , , 

^(|— ,-^)(-8y -7— +Sy ^ + i ) 
„ 2 ' 2 „ 2 Z o X ’o X 

Z Z o 

o o 
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2 


3 


3x^Z 


“1 7 4 i 

~(-4 + 20y ^-10— +4-S-5-—+20— - 1 6y 
4 o X 3 2 X o 2 


+ -i-j-(32y -96y ’ + 82— ---17-| +-i.24 -2, 

42^ ^ ° * * x’ x’ x’ 

o 

3 4 

y y y 

+ 49 "T” - 16-^ - 104 — ^ + 64y ^ ) 

2 4 X o 

X X 


with the initial conditions 


Z (x ) - 1 ) Z.(x ) — Z-(x ) ~ > • • 

O O 1 O 4 0 


The integrations of Eqs. (53) is accomplished by successive 
quadratures and its success depends on whether or not the integrals can 
be expressed in terms of known functions. It has been found that the 
following recurrence formula is very useful. We have 

J p(x)y dx = - ^i^y ” + / p(x)y ”’^dx + J ”dx 

*' o n o •'o Lx n J o 


where n 0 and p(x) is any arbitrary function. To derive the formula 
we use the well-known relation 


xl (x) + nl (-X) = xl (x) 

n n n- 1 


so that for n = 1 




and for n = 0 


I ' = I,(x) . 

o 1 


Therefore, if y =1/1, 
o o I 
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1*11/ y 

, _ Q O 1 _ j ^ ‘ 

x'^o 

* ^1 


Now consider 


/p(x,d(^) = £S2i y" . f il<2Ly“dx 

n o ^ n 'o 


/ ^0”'^ V *** ■ / p^*^ ^0”'^ ^ ‘ 


X ' o 


POOy " 

n o 


/ £lWy ■> dX . 
n o 


Rearranging the eqtiation, we have the recurrence formula (55). 

Using these relations, we proceed with the integrations of the 
Eq. (53) using the initial conditions (54). 


We have 


and by Eq. (57) 


(X) » 1 

o 


Xlj (X) 

3t_I,(x ) 
O 1 o 


where x = 3a e is the initial value of x. 
o ^00 

With Z =1, the equation for 2, is 

O Z 


2x + y - 2xy 
o o 


Integrating 


Z, = x“ + 


r 2 

Log xlj(x) - Z j • 


By the recurrence formula (55) with p(x) = x, n = 1 


f xy '“dx - -tx** - xy + 2Lo«xI, (x) 

^ ' n 2. a ■'1 
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•o that, in £q. (62) with the initial conditions, we have 

xl (x) 

Zj • 2xy^(x) - 2x^y^(x^) . 3 Log ^ ^ y • (64) 

o 1 o 

The integrations for obtaining Z^. and are performed the 
same way, but they are much more laborious. It is found that the 
Z.(x) can be expressed in terms of two functions 

I (x) X I (x) 

We have the final solution 


Z (X) 

o 

Zj(x) 


= 1 


= Log 


xlj(x) 
o 1 o 


z^Cx) 


ZjCx) 


Z4(x) 


= 2(A-A^) - 3Zj 

= j - ^(A-A^) - 2(A^-A^^) + 1 3Zj - 2AZj + 1 Zj^ 

35 22 71 22833 

= - ~(x^.x + -^(A-A J + 3{A^-A^ ) + 7 (A -A ) 

2 o 2 O 0^0 

+ 4A (A -A ) - 2(x^A-x ^A ) 

o o o o 


Z5(x) 


i^z ^ 
2 1 


-(69+ 6A +7x^ - 19A - 4A^Z, 

O 1 

-Zj^ + 2AZj^ 

2 41 3 3 4 

= Z, (162+6A ) + ^Z, Z, 

1 0 4 14 1 


+ Z,(437-~x +6A - 2Z, ^A - 6Z,^A^ 


'1 

69 


2 o 


1 


1 


- 2 Zj^A+-^Zj^x^- 8A^Zj.21A^Zj +6x^AZ^ 

X „ y -343 ^ , 147 2^ 

+AZ ( 8A ) + —X Z 

14 0 4 1 
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( 66 ) 


+ t(x^-x ^ + (14A +^)(x^.x ^ 

+ (7x/ - 39 A . 4A„^ . (A-A„) 

O 0 0 4 o 

+ 4x^A^- 4x ^A ^ + 2A^-2A 4A^ + 4A ^ . 

o o o o 

The 8emi>major axie of the orbit under contraction i* 

~ = 1 + € Zj + +€ . (6 

Using X as a parameter, we easily express the other quantities 
of interest. The eccentricity e is given by Eq. (49) while the drop in the 
periapsis is obtained from 


r -r 

p r 1 

H — * Pj^p -a(l-e)J * pa^- pa^e^+pae - p a^(l+€ Zj+. ..) 


r -r 

Po P 


2 3 4 

= (x-x^) . (Zj +€ Z^fC Z^+€ Zg). 


The ratio of the eccentricity can be obtained from 


— = - (— ) . 

e Z ' X ' 

o o 


For each initial value € = H/ a , and eccentricitv e , we can 

o ' o 

calculate the initial value x = & a e - e /€ and evaluate the expressions 

o o o o 


fora/a , e/ e and (r -r )/ H ae a function of x. Then these functions 
o o P P 

i-p 

can be cross plotted in any combination. 


The orbital period is 



( 70 ) 


L a ^ , (Z(X))^^^ . 

■*■0 o 

Since the qxientitiee a, e, and T are all easily observable, and the 

4 

integration has been performed to e included, for small and moderate 
eccentricity, which is the case for most scientific Earth satellites, the 
equations can be used to verify the assumption made on the atmosphere. 

In general, the density, as a first approximation, can be assumed to be 
locally exponential. That is to say, the parameter |3 , or H s 1/ ^ , 
can be assumed constant for each layer of the atmosphere. Since the 
value of j3 enters the analytic solution, by adjusting for concordance 
between theory and observation, p can be determined. 

The expansions used to obtain the basic nonlinear equation (51) 
are only valid for small and moderate values of eccentricity. To be 
exact, expansions in elliptic motion apply to eccentricities which are 
less than 0.663. Above that value, the series are no longer absolutely 
convergent. King-Hele used the equation, truncated to the order of e^ 
and showed that his theory is accurate for orbits with e less than 0. 2. 
Because of the difficulty he encountered in the analytic integration, he 
divided the contraction of the orbit into two phases. For phase 1 , the 
eccentricity has the approximate range 0.02 < e <0.2. In this phase 
the range of x is 3<x<30 and the integration is performed using the 
asymptotic form of the Bessel ratios y^(x). For phase 2, the eccentricity 
has the approximate range of 0 < e < 0.02, and the range for x is 
0 < X < 3. Simplification by asymptotic expansion is not available, but 
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integration to the first order in e is feasible. Also, for the case of 
large values of x, the integration by King-Hele involves a heuristic step 
in that he assumes a certain form for the solution Z on the right hand 
side of £q. (51) in its asymptotic form so that finally Z can be obtained 
by quadrature. 

For the present formulation, the nonlinear Equation (51), to the 

4 5 

order of e has been integrated rigorously to the order of € . The 

artificial division of the contraction of the orbit into two phases has been 

removed and the solution is uniformly valid for the entire lifetime of 

the orbit for e in the range 0 < e < e^ . 

Explicit Formulas For The Orbital Elements 

For small and moderate eccentricity (e < .4). the solutions obtained, 

Eqs. (66)-(70), are very accurate. This has been verified by computing 

the numerical integration of the nonlinear equation (51) and its analytic 

solution, Eqs. (66) and (67). In this case, the solution was always fovind 

to be greate ^ numerical integration with a maximum error of 

approximately ce 5(l-e ^) for 0. 1 <e < 0.99. It is interesting to 
o o o • 

note that even as e^ -• 1 , which is outside the region cf strict mathematical 

validity, the maximum error is less than 1 / 10 p r (of the order 10 ^). 

^o 

For small values of eccentricity the solution is extremely accurate. For 

example, when e =0.1 and e = 0.008. Eqs. (o6) and (t>7) provide 7 digits 
o 

of accuracy, while for the same case King-Hele' s solution gives only 4 
digits of accuracy. Thus, the present solution provides a major improve- 
ment in that it is much more accurate and that it is uniformly valid for 
all eccentricities, 0 < e < e . 



The solutions obtained for the orbital elements, a/ a , e/ e , 

o o 

etc. are expressed in terms of the variable x. It would be useful to derive 
formulas to relate any pair of orbital elements. This amounts to the 
elimination of x between any two of the equations (67) - (70). Because 
of the transcendental nature of the solutions, the task is cumbersome. 
Fortunately, since € is small, elimination of x through Lagrange* s 
expansion is feasible. 

To derive the explicit expression for the semi-major axis in 
terms of the eccentricity, we write Eq. (67) as 


Z = p + c (j) (Z) 


(71) 


where, by observing that 


X = ^ a e = 


we can write 


X = o Z 


with 


8 a e (^) ( ~ ) 
o o a e 
o o 


a - X (— ) 
o e 

o 


(72) 


(73) 


Then, explicitly 


b (Z) = Z.(a Z) + € Z (a Z) + c Z (a Z) + ... 

i ^ 3 

If Lagrange* s expansion is applied to Eq. (71), we have 

00 n n-1 

Z = p+ E ^(;~) [6(p)j". 


n=l 


nl dp 


(74) 


(7 5) 


If we carry out the expansion, and then put p s 1 , we shall have to the order 

of e® 
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(76) 


Z s 1 + € hj( a ) + € **^2^**^ ^ ^ '> ■*■ f ^ + € ^^5^® ) 


where 


h. s 


h 

h- 


= 2(A-A ) + (A-3)Z, 

o i 

= lia^.x ^) - (2A + ~)(A..A )+(l3+2a^.4A-A^)Z, 

Z Q 0 4-0 1 

12 2 2 

+ -^{3+0 +A-A )Z^ 

= - 4(a^-x ^)(35+4A -7A) 

4 0 o 

+ i(A-A H213+42A +16A ^+12o^-9A+4A A-8A^) 
bo o o o 

- J Zj(138+25a^-46A-7A^-2A^) 

- ~ Z ’(35+7o^-6A^+o^A-A^) 

im X 

+ 2(A-A )Z,(3+ff ^+A-A^) 

o r 

- i Z,^(6-o^+2a^A + 3A^-2A^) 

t> 1 

322^122 2 ’ 

= - y(« -X -X )(885+b8o -92A-12A“) 

4 O o O 

-(o^-x ^HA-A H-^+4A +2A) 

O 0 4 0 

,, , , 441 33 2 137 ^ , 2, 9 ^2 ,3, 

-(A-A )( -V+ -Tro + — r— A+2a A + — A +2A ) 
o 2 2 4 2 

2 409 2 2 3 40 

+(A-A ) {^+2a +23A+12A )-(A-A ) (104 — A) 
o 8 o 3 


(77) 


44(A-A ) 4 ■r(a -X )Z (3+a +A»A ) 

o 2 o 1 

+ 2 ^( 437 + 880 ^+ 20 **- 1 54A-5o ^A- 22A“-3A^-A**) 

- t(«-A )Z,(143+29o"+25A + 12o“A.17A“-12A^) 

2 o I 

+ 2(A-A )“Z,(3+o^+A-A”) 
o 1 

+ iz,“(b48+133o“+4o**-96A+14o“A-91A“-4a 'A"-1:A i 
4 I 

- 'A-A )Z."(b.o"+2o‘'A + 3A"-2A^) 

o 1 

13 ^4 ’’ ’22 34 

4 - Z, (12344o‘'-o -6A4l3o A + 15A 42o A -ISA -A ) 

b 1 

14 74 .? 777 34 

4— Z, (18-o"-2o -80 A-3A"+8o‘A 412A -oA ) 

24 1 


-2b- 


o> 


In the expressions for h^(a), we have defined 

al^(a) 

A . a = Log 

o 1 o 


( 78 ) 


Since or = x (e/ e ), the solution, as given by Eqs. (76) and (77) provides 
o o 

an explicit expression of the variation of a/ a as a function of e/ e . 

o o 

The other orbital elements can also be expressed in terms of e. 
The drop in the periapsis is 


r - r 
P P 

— = (a-x^) - (1-e) (hj+ch^+c ^h^+€ ^h^+€^hg) . 


(79) 


For the apoapsis, we have 


r -r 
a a 
o 


H 


= (x-o) - (l+e)(h +€h +e^h,+c^h +c^h ) 

o 14 3 4 5 


(80) 


In addition, the following formulas can be easily derived 


r 

P 


r 


P 


o 


(1-Ca) 

(1 - .„) 


Z(«) 


r 


a 


r 

a 

o 


(H-eg) 

(l+e ) 
o 


Z(ai 


(81) 


(82) 


_T 

T 



(83) 


o 

The last expression provides the orbital period as a function of the 
eccentricity. As pointed out by King-Hele, such a relationship, if it is 
accurate, which is the case in the present theory, provides a powerful 
method of verifying the assumption made on the atmosphere from two of 
the most accurate and easily measured orbital parameters, namely, the 
period of revolution and the eccentricity. 
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The Contraction of tiighly Eccentric Orbits 


For the case of orbits with large eccentricities, King-Hele has 
used an entirely different method to derive the basic equation for the 
contraction. Using the present notation, we have his eqxiation in the 
form 


H = ■ ZTTT 


(84) 


It is possible to obtain this equation from the basic equation (51). 

Since x = pae, when e -• 1, a -• oo , x becomes very large and the asymptotic 

expression for Bessel' s ratio y (x) is 

o 


+ ••• 

y ^ = 1 + i + ... 

O X 


(85) 


Using this form in Eq. (51) we have 


dZ 1 


3 4 2 5 3 

€ X € X ex 


( 86 ) 


an equation whivji? can be seen as the development of Eq. (84). While 
King-Hele can only provide an approximate solution to the nonlinear 
equation (84) by assviming that on the right hand side Z is approximated 
by Z = 1 - e the equation can be integrated by simple 

quadrature, it is enlightening to know that the exact solution to the 
equation can be obtained. This is done by using the transformation 

z = e (X +q) (87) 

and changing the independent variable from x to q to have the Bernoulli 
equation 
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dq 


( 88 ) 


= 2x + 


4x 


Using the change of variable, 

_ 2q 


X = 


K(q) 


( 89 ) 


we have after substituting into (88) 

2q 


dK 

dq 


-4e 


( 90 ) 


Integrating, we find that K can be expressed in terms of the exponential 
integral. 


K = -4/ 


2q 


2q 


•d (2q) + C 


*4 C*(2q) 4* C . 


Thus, we have the exact solution in parametric form. 


2(q-qQ) 


X = 


-2q 


^ + 4e ° [Ei(2q^> - E.(2q)] 


along with Eq. (87) and the initial conditions 


( 91 ) 


( 92 ) 


Z(Xo) = 1 


= — - X 

€ o 


(l.eo) 


(93) 


To evaluate the exponential integral, we first observe that the 
argument 2q is large, hence its asymptotic form is adequate. In general, 
we consider the integral 

En(x) ■ . (94) 

By integration by part 

« ^ . n-1 X 

E^(x) = X e - (n-1) E^_^ 



By repeated application ol this formula, we deduce the asymptotic expansion 


for large x 


E (X) = [ 1 . te-.i)(»-Z)(n-3) ^ ... ] 

XX ■ 

'When n=0, we have the exponential integral and by taking 6 terms of the 


(95) 


series, for x > 5, the solution is identical to the numerical values tabulated 
by Abramowitz and Stegiin (1972). 


Numerical Application 

The accuracy of the present theory has been verified by comparing 

the analytic solution and the exact numerical solution for a wide range of 

orbits about the Earth' s atmosphere. The numerical solution is obtained 

by integrating the basic nonlinear equation (51). The analytic solution 

employed is either the solution Z = Z(x) as given by Eqs. ( 66 ) and (67) or 

the explicit solution Z = Z(a) as given by Eqs. (76) and (77). 

The parameters used are the initial eccentricity e^ and the initial 

perigee distance r , or equivalently the dimensionless small parameter 
^o 

1/ 6 r . Then we have 
P 

( 1 -eJ 

€ = 




(96) 


which tends toward zero as e 1 . By using the three values — 

o / o 0 j 


Sr 

Po 

= 0.005, 0.01 and 0.02 we cover a wide range of perigee heii’hts. 

Figure 3 plots the variation of Z = a/ a^ as function of x/x^ for 
different values of the eccentricity. Since the analytical solution has a 
high degree of accuracy, its small deviation in the fifth or sixth digit 


from the nvimerical solution cannot be detected in the figure. Initially 

the solut'.'n is nearly linear with the slope in the figure approximately 

equal to e^. but near the end. as x and e approach zero, it exhibits rapid 

decay. This explains the difficulty encountered by King-Hele in his analytic 

integration. It is interesting to note that the analytic solution Z(x) of the 

trxuicatcd eq\iation (51) remains accurate for high eccentricity. For e = .99, 

o 

this solution and the exact solution (87) and (92) of the asymptotic equation 

(84) are nearly identical except for very small values of x / x^ . 

Figure 4 plots the solution Z = Z(ar) = Z(e/ c ) as a function of the 

eccentricity for several values of e^. The range of validity is limited to 

e^ = 0. 5 since the Z(ar ) solution is not as accurate as the Z(x) solution. 

It was found that Z(or) exceeds the ntimerical solution by a maximum value 

approximated by e 15(l-e ), which still gives 7 digits of accuracy for 
o o 

e^ = 0. 1 , but diverges as 1* For e^ = 0. 5 the error is imperceptible 

in the figure. The decay in the perigee distance r^ versus the eccentricity 

for e = 0. 1 and e = 0.4 is presented in Fig. 5. The ratio r / r remains 
o o ® P P 

o 

nearly equal to one for a large portion of the decay process, but as e -► 0 

the drop in perigee altitude increases rapidly. For small values of 
1 


p r 


(0. 005) the decay is much slower than for large values (0. 02) as 


Po 


evident in the figure. The fractional error in this plot and in the next two 
is kept to an imperceptible amount by considering the error formula 
mentioned earlier. 

The decay in the apogee distance r versus the eccentricity for 

& 

several values of e^ is presented in Fig. 6. It is evident that the ratio 
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r / r decreases rapidly with the eccentricity. Initially the parameter 

A el 

o 

■■ ' seems to have little effect, but as the eccentricity approaches 

aero the larger values of yield more rapid decay. 

Finally the decay in the orbital period T as function of the eccentricity 
is presented in Fig. 7. As pointed out by King-Hele this functional relation- 
ship T = T{e,€) provides a powerful formula for testing the atmospheric 
parameter c since orbital period and eccentricity can be accurately 
measured. 
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Fig. 1. Nomenclature. 
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Fig. Z. .Aerodynamic Forces . 
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